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ABSTRACT 

The solar chromosphere contains thin, highly dynamic strands of plasma known as spicules. Recently, it has 
been suggested that the smallest and fastest (Type II) spicules are identical to intermittent jets observed by the 
Interface Region Imaging Spectrograph. These jets appear to expand out along open magnetic field lines rooted 
in unipolar network regions of coronal holes. In this paper we revisit a thirty-year-old idea that spicules may 
be caused by upward forces associated with Alfven waves. These forces involve the conversion of transverse 
Alfven waves into compressive acoustic-like waves that steepen into shocks. The repeated buffeting due to 
upward shock propagation causes nonthermal expansion of the chromosphere and a transient levitation of the 
transition region. Some older models of wave-driven spicules assumed sinusoidal wave inputs, but the solar 
atmosphere is highly turbulent and stochastic. Thus, we model this process using the output of a time-dependent 
simulation of reduced magnetohydrodynamic turbulence. The resulting mode-converted compressive waves are 
strongly variable in time, with a higher transition region occurring when the amplitudes are large and a lower 
transition region when the amplitudes are small. In this picture, the transition region bobs up and down by 
several Mm on timescales less than a minute. These motions produce narrow, intermittent extensions of the 
chromosphere that have similar properties as the observed jets and Type II spicules. 
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1. INTRODUCTION 

The Sun’s hot corona expands into interplanetary space as 
a supersonic plasma outflow known as the solar wind. How¬ 
ever, we still do not yet know how the tenuous corona/wind 
system is formed from the much larger pool of mass and 
energy in the colder photosphere and chromosphere. High- 
resolution observations (e.g., Eletcher et al. 2015) show that 
coronal heating is highly dynamic and intermittent in space 
and time. Some have suggested that much of the corona’s 
mass and energy may be injected in the form of narrow fea¬ 
tures known variously as spicules, fibrils, and mottles (Beck¬ 
ers 1972; Sterling 2000; De Pontieu et al. 2007; Tsiropoula 
et al. 2012). Recently, the Interface Region Imaging Spectro¬ 
graph {IRIS, De Pontieu et al. 2014a) found similar jet-like 
features emerging from largely unipolar network flux concen¬ 
trations in coronal holes and quiet regions on the solar surface 
(Tian et al. 2014, 2015). 

Coronal holes have long been known to play host to nar¬ 
row, ray-like features known as polar plumes and polar jets 
(Newkirk & Harvey 1968; Ahmad & Withbroe 1977; Wang et 
al. 1998; Dobrzycka et al. 2002; Culhane et al. 2007; Raouafl 
et al. 2008; Young & Muglach 2014; Paraschiv et al. 2015). 
These bright strands typically extend up to heights of order 
0.1-1 Rq, whereas spicules and the IRIS network jets have 
length scales of only ~0.01 Rq. The largest plumes and 
jets are often associated with the emergence of small mag¬ 
netic loops at their footpoints, and thus they are believed to be 
powered by magnetic reconnection (see, e.g., Wang & Shee- 
ley 1995; Shibata et al. 2007; Pariat et al. 2009; Moore et al. 
2010; Cheung et al. 2015). However, no clear evidence for 
reconnection in the smaller IRIS jets has been found so far. 
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In this paper we explore the idea that magnetohydrody¬ 
namic (MHD) waves are responsible for producing rapidly 
varying, field-aligned extensions of cool chromospheric gas 
along network flux tubes in coronal holes. These modeled 
features are found to have similar properties as the IRIS net¬ 
work jets, which Tian et al. (2014) proposed to be identical to 
the so-called Type II spicules observed above the limb. The 
idea that spicules could be driven by MHD waves has been 
discussed for several decades (Hollweg et al. 1982; Mariska 
& Hollweg 1985; Hollweg 1992; Kudoh & Shibata 1999; De 
Pontieu 1999; Velli & Liewer 1999; Matsumoto & Shibata 
2010; Murawski et al. 2015), but it is still not known whether 
this is a dominant mechanism in the real solar atmosphere. 

Many of the necessary ingredients of the wave-driven 
spicule model are supported by observations. The solar 
corona contains transverse, incompressible oscillations in the 
magnetic field and plasma velocity (Banerjee et al. 1998; 
Tomczyk et al. 2007; Jess et al. 2009; McIntosh 2012), but 
whether they should be called “Alfven waves” is still a mat¬ 
ter of debate (e.g., Mathioudakis et al. 2013). More specifi¬ 
cally, spicules and jets themselves seem to contain torsional or 
kink motions that could have an Alfvenic character (Kukhian- 
idze et al. 2006; De Pontieu et al. 2014b; Tavabi et al. 2015). 
There is also evidence for longitudinal compressive waves— 
i.e., fluctuations in density and the velocity component par¬ 
allel to the wavenumber vector—that may or may not follow 
the ideal MHD magnetoacoustic dispersion relations (Ofman 
et al. 1999; Krishna Prasad et al. 2012; Threlfall et al. 2013; 
Miyamoto et al. 2014). In polar plumes, Liu et al. (2015) 
found high-frequency Alfven-like waves and low-frequency 
compressive waves traveling along the same held lines, which 
points to the possibility of mode coupling. Also, Pant et 
al. (2015) found what appear to be slow-mode magnetosonic 
waves within the IRIS network jets. 

The models developed in this paper rely on MHD waves be¬ 
having in an intermittent and stochastic manner. Several ear- 
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Her attempts to understand spicules as a by-product of waves 
assumed a periodic, sinusoidal driver at the lower boundary of 
the modeled system. However, both observations (Tomczyk 
& McIntosh 2009; Liu et al. 2014) and simulations (e.g., van 
Ballegooijen et al. 2011; Perez & Chandran 2013; Zhdankin 
et al. 2015) show that MHD fluctuations in the chromosphere 
and corona exhibit continuous power-law spectra and irregu¬ 
lar bursts of activity and dissipation. This variability may also 
be related to the fact that the Sun’s transition region (TR) has 
a complex “corrugated” shape (Feldman et al. 1979; Zhang 
et al. 1998; Peter 2013). We propose that jets and Type II 
spicules are short-lived extensions of the corrugated TR that 
are driven by similarly infrequent outliers in the underlying 
population of waves and turbulent eddies. 

The remainder of this paper is organized as follows. Sec¬ 
tion 2 discusses how chromospheric Alfven waves may evolve 
nonlinearly into a collection of compressible fluctuations. In 
Section 3 we estimate the degree of nonthermal expansion 
experienced by the upper chromosphere and TR as a result 
of shocks that develop from a compressive wave train. Sec¬ 
tion 4 takes the output from a reduced MHD simulation of 
Alfvenic turbulence in a coronal hole and computes the time- 
dependent generation of compressive waves and intermittent 
levitation of the TR. In Section 5 we compare the modeled 
up-and-down motions of the TR with the observed properties 
of Type II spicules and IRIS network jets. Lastly, in Section 6 
we summarize our results, discuss some broader implications, 
and describe future improvements to the models. 

2. NONLINEAR PRODUCTION OF PAR ALT FT. VELOCITY 
ELUCTUATIONS 

When incompressible Alfvenic fluctuations grow to a suf¬ 
ficiently large amplitude (i.e., when the oscillating transverse 
magnetic field SB± becomes of the same order of magnitude 
as the background held strength Bq), they become susceptible 
to a range of nonlinear interactions that can spawn other types 
of waves. For example, 

1. It has been known for several decades that linearly po¬ 
larized Alfven waves can excite second-order pondero- 
motive oscillations in density, gas pressure, and mag¬ 
netic pressure (Hollweg 1971; Spangler 1989; Vasquez 
& Hollweg 1996). Because these oscillations are tied 
to the extrema of transverse arcs traced by the magnetic 
held vector, their frequencies tend to be twice those of 
the original Alfven wave. Corresponding wave periods, 
for Alfven waves oscillating with Va ~ 3-5 minutes, 
are of order 1-2 minutes. These are reminiscent of the 
durations and recurrence timescales of Type II spicules. 


2. There has also been substantial work done to study 
the nonlinear development of parametric instabilities 
for circularly polarized Alfven waves (Goldstein 1978; 
Jayanti & Hollweg 1993; Turkman! & Torkelsson 2003; 
Del Zanna et al. 2015). For conditions appropriate 
to the corona, this instability usually involves an up¬ 
ward propagating Alfven wave decaying into a down¬ 
ward propagating Alfven wave and an upward propa¬ 
gating magnetosonic-like wave. The latter tends to have 
a lower frequency ws than that of the original Alfven 
wave wa- Typical periods {V = 27r/a;) in the open-field 
corona are 


Vs 


Va 


( 1 ) 


which for Va ~ 3 minutes, typical coronal Alfven 
speeds Va ~ 2000 km s“', and sound speeds ~ 150 
km s“', gives a period of order 20 minutes, similar to 
what is observed for density fluctuations above the limb 
(e.g., Ofman et al. 1999; Threlfall et al. 2013). 


The number of possible mode-coupling interactions grows 
even bigger when the MHD waves pass through a strongly in¬ 
homogeneous background medium (e.g., Heyvaerts & Priest 
1983; Lee & Roberts 1986; Nakariakov et al. 1998; Hollweg 
& Kaghashvili 2012; Thurgood & McLaughlin 2013). 

Because we are concerned with the development of short¬ 
lived spicules and jets along unipolar held lines, we focus on 
the most rapid kind of nonlinear mode coupling that occurs 
in homogeneous plasmas: the second-order ponderomotive 
effect. Hollweg (1971) showed that large-amplitude linearly 
polarized Alfven waves produce an oscillation in the par¬ 
allel velocity with an amplitude that scales as 



SBa^ 


2 


( 2 ) 


where Np describes the dependence on the plasma /3 param¬ 
eter. We use a simplified kinetic definition of /3 = (Cj/Va)^ 
which is different by a factor of 2/7 Ri 1.2 from the standard 
MHD definition of the ratio of pressures. Hollweg (1971) 
used the ideal MHD conservation equations to derive 


Np 


0.25 


(3) 


which diverges unrealistically at the asymptotic value of /3 = 
1. Vasquez & Hollweg (1996) extended the second order 
MHD theory and found, for some cases, Np = 0.25/(1-F/3). 
Several other properties of this solution were consistent with 
a magnetosonic-like wave mode. 

A substantial amount of other theoretical work has been 
done to simulate the nonlinear evolution of Alfven waves, 
much of which involves solving the Derivative Nonlinear 
Schrodinger (DNLS) equation (e.g., Mjolhus &. Wyller 1986; 
Medvedev & Diamond 1996). In this paper, we make use 
of the kinetic results of Spangler (1989), who solved a per¬ 
turbed Vlasov equation for the ponderomotive density fluctua¬ 
tions associated with a soliton-like Alfvenic pulse. We solved 
the Spangler (1989) equations for a range of plasma (3 values 
and for the simple one-temperature case of Te = Tp. Effective 
compressional amplitudes were extracted from the simulated 
density fluctuation profiles, which did not maintain the same 
Lorentzian shape of the input Alfvenic pulse. We maintained 
continuity with earlier studies of sinusoidal waves by comput¬ 
ing the relative fluctuation ratio Spjpo as half the peak-to-peak 
pulse variation in density. Also, we assumed the acoustic-like 
energy equipartition found by Hollweg (1971), 


, 4 ) 

Va Po 

and used it to compute (5v|| as in Equation (2). The numeri¬ 
cal results were used to construct the efficiency factor A(a, and 
Eigure 1 shows its dependence on f3. This function agrees 
with the above analytic results in the limits of /3 <C 1 and 
/3 ^ 1, and it falls in between them for /3 « 1. We fit the 
P dependence of this factor by an approximate function. 


0.25 , 0.135^2-4 

0.305+ /34.6 


2c, 


(5) 
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Figure 1. Dependence of the dimensionless mode-conversion efficiency 
on the kinetic plasma /3, including results from Hollweg (1971) (dot-dashed 
curve), Vasquez & Hollweg (1996) (dotted curve), and present work based on 
Spangler (1989) (solid curve). 

and this is accurate to within about 3% over the range of /3 
values shown in Figure 1 . A plot of this function would be 
nearly indistinguishable from the solid curve that shows the 
numerical results of the Spangler (1989) model. 

Later in this paper, Equations (2), (4), and (5) are used to 
compute the properties of acoustic-like waves that are gener¬ 
ated from a time-dependent simulation of Alfvenic turbulence 
in the chromosphere and low corona. We find that there is a 
distinct local maximum in (5v'|| at a height Zb = 876 km in the 
low chromosphere. That location will be taken as a concep¬ 
tual “base height” and used to compare against other models 
of acoustic wave evolution. 

3. SHOCK STEEPENING AND CHROMOSPHERIC LEVITATION 

In the solar chromosphere, upwardly propagating acoustic 

—1 /2 

waves undergo rapid growth in amplitude, with (5v|| cx Pq 
in the limit of an isothermal atmosphere and no dissipation 
(Lamb 1908, 1932). At some point, however, the amplitudes 
become large enough for the waves to steepen into shocks 
(e.g., Stein & Schwartz 1972) and thereby dissipate their en¬ 
ergy into the surrounding atmosphere. In the process of damp¬ 
ing, the fluctuations are also able to exert a mean upward 
wave pressure gradient force on the atmosphere (Dewar 1970; 
Jacques 1977). This force is essentially a net transfer of mo¬ 
mentum from the oscillations to the background gas. In a 
hydrostatic atmosphere (or in the subsonic parts of the so¬ 
lar wind), an increase in the total effective pressure increases 
the gravitational scale height, and this in turn “puffs up” the 
cool chromosphere. This effect has also been explored in the 
context of atmospheres of cool, evolved giant stars (see, e.g., 
Bertschinger & Chevalier 1985; Ludwig & Kucinskas 2012). 

We simulated the above chain of events using a series of 
one-dimensional time-steady solutions of the ZEPHYR code 
(Cranmer et al. 2007). This code produces a realistic de¬ 
scription of the photosphere, chromosphere, corona, and solar 
wind in open-field regions of the solar atmosphere. ZEPHYR 
solves equations of wave action conservation for both acous- 
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Figure 2. Height dependence of (a) time-steady temperatures T and hydro¬ 
gen number densities nn, and (b) acoustic wave velocity amplitudes (5v||, 
corresponding to a set of ZEPHYR models with a range of acoustic wave 
fluxes Fs. Line colors/styles denote Fs (see caption) and are consistent in 
both panels. Also shown are the chromospheric base height Zb (white circles) 
and transition region height ztr (black circles). 

tic and Alfven waves in the presence of several kinds of dis¬ 
sipation and shock steepening, and it includes wave pressure 
terms that couple the fluctuations to the background plasma. 
In these models there is no coupling between the acoustic and 
Alfven wave modes; they each evolve independently of one 
another. Specifically, we made use of the grid of models from 
Section 8.2 of Cranmer et al. (2007), in which the Alfven 
wave properties were held fixed and the acoustic wave flux 
at the photospheric boundary was varied over several orders 
of magnitude. 

Eigure 2 summarizes the results from the ZEPHYR models 
with a range of acoustic wave power inputs. Eigure 2(a) shows 
the time-steady variation of density and temperature with 
height, and Eigure 2(b) shows how the root-mean-squared 
(rms) parallel velocity amplitude varies with height for 
these models. As the acoustic waves become stronger, the 
chromospheric scale height receives an increasingly large 
augmentation from wave pressure. In the models of Cran¬ 
mer et al. (2007), the sharp TR between chromospheric and 
coronal temperatures occurs when the density dips below a 
critical value determined by the peak of the optically thin ra¬ 
diative loss function (see also Owocki 2004). The models 
with stronger acoustic waves have flatter density gradients, so 
the critical density is reached at larger heights. 
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Table 1 

Chromospheric Properties of ZEPHYR Models 


Fs 

(erg cm“^) 

<5v|| (Zi) 
(km s“') 

<5v|| (ztr) 
(km ) 

ZTR 

(km) 

UlAU 
(km ) 

M 

(MqIjt) 

0 

0 

0 

3364 

723 

1.91 X 10-“* 

10^ 

0.0416 

4.22 

3508 

722 

1.91 X lO^*'* 

10'’ 

0.132 

1.11 

4078 

724 

1.88 X 10-*'* 

10’ 

0.416 

11.4 

5062 

720 

1.90 X 10-‘4 

10* 

1.32 

14.3 

6540 

720 

1.88 X 10-*'* 

lO"* 

3.26 

16.7 

8484 

721 

1.86 X 10“*^ 

10‘° 

5.36 

18.1 

10576 

728 

1.76 X 10-*'* 


Table 1 lists some key properties of each ZEPHYR model. 
The acoustic wave energy flux is injected at the photo- 
spheric lower boundary, and the velocity amplitudes (5v|| in¬ 
crease monotonically with increasing height from the photo¬ 
sphere to the TR. Between the photosphere (z = 0) and the 
TR height (ztr) we focus on the lower chromospheric base 
height Zb = 876 km and highlight it with a vertical line in Fig¬ 
ure 2(b). The TR is defined as the height at which the modeled 
temperature T first rises to 20,000 K. This is lower than the 
temperatures at which most TR emission lines are formed, but 
we are using ztr as an effective height for the tips of spicules 
and jets; i.e., where chromospheric emission ends. 

The time-averaged shock-driven levitation of the TR is a 
key result of the models shown in Figure 2. Without any com¬ 
pressive waves, the Sun’s TR occurs in the models at a height 
of ^^3300 km above the photosphere. This is slightly higher 
than the canonical range of 2000-2500 km that is seen in 
one-dimensional empirical models (e.g., Vernazza et al. 1981; 
Avrett & Loeser 2008). With increasing “turbulent pressure,” 
the modeled ztr can increase to values greater than 10,000 
km. Table 1 gives a span of heights between these extremes, 
and these are roughly consistent with the values sometimes 
reported from off-limb chromospheric measurements (Zhang 
et al. 1998; Filippov & Koutchmy 2000). In any case, we use 
the results given in Table 1 as an interpolation lookup table 
that provides an instantaneous estimate of ztr for any given 
value of (5v||(zfe). 

In the upper chromosphere, the acoustic waves begin to dis¬ 
sipate because they steepen into shocks. The time-averaged 
damping rate is dominated by the entropy change T AS at each 
shock, spread out over the time between successive shock pas¬ 
sages. The models with larger values of Fs undergo steepen¬ 
ing at lower heights, so their chromospheric profiles of i5v||(z) 
become flatter and more saturated. Above the TR, the acoustic 
waves damp rapidly due to the greatly amplified rate of heat 
conduction in the corona (see, e.g.. Equation (26) of Cran- 
mer et al. 2007). Because of this rapid damping, the acoustic 
waves have a negligible effect on the eventual acceleration of 
the fast solar wind.^ Note from Table 1 that the wind speed 
at 1 AU and the sphere-averaged mass loss rate M are barely 
affected by changing Fs. 

Below, we make use of the correlation seen in Table 1 be¬ 
tween (5v'[|(zi) and zxr- This is applied to a time-dependent 
simulation of incompressible turbulence that does not con¬ 
tain acoustic waves driven at the photosphere. Instead, we 

^ However, Nuevo et al. (2013) suggested that the rapid damping of com- 
pressive waves may have a strong impact on the coronal heating in large hy¬ 
drostatic loops. With enough energy dissipated at low heights, the maximum 
temperature may occur below the loop apex, thus giving rise to loops with 
decreasing T{z) in the corona. 


assume that compressive waves are generated throughout the 
chromosphere by the nonlinear mode conversion mechanism 
discussed in Section 2. We assume that the upward evolution 
of those compressive waves produces a shock-driven levita¬ 
tion similar to that seen in the ZEPHYR models. Of course, it 
should be made clear that the second-order variations in (5v|| 
described by Equation (2) are not identical to classical sound 
waves; e.g., they propagate at a phase speed of Va instead 
of Cj. However, once these waves steepen into shocks, the 
time-averaged loss of momentum and energy is expected to 
be similar to the acoustic-wave case. 

4. RESULTS EROM TIME-DEPENDENT MHD TURBULENCE 

Although the ZEPHYR code simulates the transport, cas¬ 
cade, and dissipation of Alfvenic turbulence, it does so us¬ 
ing time-averaged phenomenological equations. These equa¬ 
tions do not self-consistently simulate the actual process of 
an MHD cascade, which is believed to be a consequence of 
partial wave reflections and nonlinear interactions between 
Alfven wave packets. In order to more accurately model 
these processes, a time-dependent and three-dimensional ap¬ 
proach is needed. We used a reduced MHD (RMHD) code 
called BRAID (van Ballegooijen et al. 2011; Asgari-Targhi 
& van Ballegooijen 2012; Asgari-Targhi et al. 2013; van Bal¬ 
legooijen et al. 2014) that simulates the generation and evo¬ 
lution of incompressible turbulence along an expanding flux 
tube with a circular cross section. This code has success¬ 
fully simulated the intermittent and dynamic heating seen 
in the chromospheric and coronal regions of closed loops. 
Also, Skogsrud et al. (2014) proposed that this type of tur¬ 
bulence model may explain the complex multi-threaded dy¬ 
namics seen within Type II spicules. 

We used the simulation of an open coronal-hole flux tube 
developed by Woolsey & Cranmer (2015). The BRAID coro¬ 
nal hole model extends from the solar photosphere (z = 0) 
to a maximum height of z = 2Rq. The choice of the latter 
value was a compromise between wanting to model as much 
of the solar wind’s acceleration region as possible and the fact 
that the RMHD equations in BRAID do not yet include the 
background outflow speed (i.e., they assume the wind speed 
u ^ Va. which breaks down above a few solar radii). The 
model was run for 2300 s of simulation time, which cor¬ 
responds to about three times the Alfven wave travel time 
through the radial grid. 

The background properties of the BRAID model (e.g., Bq, 
Va, Cs) are the same as in the time-steady polar coronal hole 
model of Cranmer et al. (2007). Figure 3 shows a selection of 
these quantities from the upper photosphere (z = 100 km) to 
the low corona. Note that this is specifically a model of a ver¬ 
tical flux tube rooted in the bright supergranular network; the 
field strength Bq remains greater than 100 G until one reaches 
a height of about 1100 km. The magnetized network chro¬ 
mosphere (above the so-called “merging height” where indi¬ 
vidual intergranular flux tubes join together) can be seen in 
Figure 3 as a region of rapidly decreasing j3 with increasing 
height. 

Figures 3(a) and 3(b) show the time-averaged amplitudes of 
transverse fluctuations in velocity and magnetic field, respec¬ 
tively. As described further by Woolsey & Cranmer (2015), 
the rms averaging in these cases was performed in two steps; 
(1) the variance was taken over all 92 discrete RMHD spectral 
wave modes at each height and time step, then (2) the variance 
over simulation time was computed at each height, excluding 
the earliest time steps during which the base-driven waves had 
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Figure 3. Time-averaged plasma properties in the BRAID model of coronal 
hole network: (a) Alfven speed Va (gold dotted curve), rms tranvserve ve¬ 
locity amplitude (red dashed curve), rms longitudinal compressive wave 
amplitude (solid black curve), (b) background magnetic field strength 
Bq (green dotted curve), rms transverse magnetic field amplitude SB± (blue 
dashed curve), plasma 0 ratio (black dot-dashed curve). As in Figure 2(b), 
the base height Zb is highlighted by a vertical line. 


not yet traversed the grid. Note that the BRAID fluctuations 
do not obey ideal MHD energy equipartition; it is possible for 
the energy densities in the kinetic and magnetic fluctuations 
to be unequal to one another. It is clear that the conditions in 
the upper chromosphere (around Zb) appear to be optimal for 
strong nonlinear mode conversion; not only is ^ ~ 1, but also 

To focus on the region(s) of the chromosphere and corona in 
which the mode conversion is strongest, we applied Equation 
(2) to the time-averaged rms properties shown in Figure 3. 
The resulting height dependence of the rms (5v'|| amplitude is 
also plotted in Figure 3(a). There is a clear peak in the low 
chromosphere at a height of 876 km, with a maximum am¬ 
plitude of ~1 km s“'. This location is defined as the base 
height Zb, and we consider it as an effectively localized chro¬ 
mospheric “source” of compressive waves. 

Next we examine the highly variable and intermittent na¬ 
ture of the BRAID turbulence at a fixed height. Specifically, 
at the base height zt, the fluctuations were averaged over all 92 
spectral wave modes (i.e., integrated over the k± power spec- 
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Figure 4. BRAID-model time dependence of: (a) the simulated compressive 
wave amplitude 5v|| at the base height Zb^ and (b) the instantaneous TR height 
ZTR corresponding to the variable 5vj|. Individual maxima and minima in ;:tr 
are highlighted with red circles. 


trum), but not over time. At this height, the mean value of 
the RMHD Alfven velocity amplitude is = 5.03 km s“*, 
with a standard deviation of 2.92 km s“' and a long tail in 
the distribution that extends up to a maximum value of 13.9 
km s“*. For comparison, the Alfven speed Va and sound speed 
Cs at this height are 9.64 km s“' and 8.35 km s“', respectively. 
A similarly processed time series of the ratio SB± /Bq has a 
mean of 0.45 at this height, a standard deviation of 0.21, and 
a maximum value of 1.08. 

Figure 4(a) shows the result of applying Equation (2) to es¬ 
timate the time dependence of i5v|| at Zh- Note that this wildly 
fluctuating quantity is an actual amplitude and not the full 
time-dependence of the parallel velocity. If the input Alfven 
wave had been a monochromatic sinusoidal oscillation, the 
estimated 5v|| amplitude would have been an unchanging con¬ 
stant. Also, because of the nonlinear nature of the mode con¬ 
version, i5v|i(f) ends up having a greater relative variability 
than either ov_L(f) or 5B±^{t). For these incompressible trans¬ 
verse amplitudes, the ratios of their standard deviations to 
their mean values are about 0.5-0.6. For the computed time 
series of (5v||, this ratio is about 1.0. 
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Figure 4(b) shows the time dependence of ztr, which 
we computed via straightforward interpolation from the 
ZEPHYR model properties given in Table 1. The TR height 
varies up and down with swings of order 2 to 5 Mm, which 
overlaps with the observed range of IRIS network jet lengths 
(Tian et al. 2014). Implicit in Figure 4(b) is the assumption 
that the upper chromosphere’s response to variability in (5v|| 
is more or less instantaneous. A more accurate model would 
have to include a finite relaxation time for the wave-pressure 
levitation to take effect. Because the compressive waves travel 
at a phase speed Va, we anticipate that this relaxation time 
should be given roughly by the Alfven-wave travel time from 
Zb to ztr- This travel time is about 50 s, which is similar in 
magnitude to the recurrence time between the modeled os¬ 
cillations in ztr- Thus, even though future simulations are 
needed to verify these effects (see Section 6), we do not be¬ 
lieve they will differ greatly from the simpler estimates made 
here. 

5. INTERPRETATION: SPICULES AND JETS? 

In the ZEPHYR models discussed in Section 3, the “steep¬ 
ened” acoustic wave amplitudes (5v|| at the TR were of order 
5 to 20 km s“'. These velocities are small when compared 
to the observed velocities of Type II spicules and IRIS net¬ 
work jets. However, we see much larger apparent velocities 
when examining the upward and downward variations of ztr 
in Eigure 4(b). We compute these velocities by recording the 
minima and maxima in Ztr(0 and taking sequential finite dif¬ 
ferences between them in height (Az) and in time (Af). These 
minima and maxima are shown as red symbols in Figure 4(b). 
The apparent velocity for each jet-like event is then computed 
as Vjet = Az/Af. From the oscillatory nature of Ztr(0i it is 
apparent that there will be roughly equal numbers of positive 
(low to high) and negative (high to low) values of Vjet. We 
eventually plan to simulate /R/5-like emission-line images of 
these evolving features, but doing so is beyond the scope of 
this paper. For now, we focus on the upward motions in ztr 
(e.g., the 46 out of 93 cataloged events for which Vjet > 0) and 
compare them to observed upward motions in Type II spicules 
and network jets. 

Figure 5 shows how the modeled collection of positive Vjet 
values is correlated with their “lifetimes” Af. We found in 
the list of 46 events that Az is roughly proportional to Af^, 
so there is a roughly linear relationship between Vjet and Af. 
In Figure 5 we also show reported ranges for the speeds and 
lifetimes of Type I and II spicules (Pereira et al. 2012) and 
IRIS network jets (Tian et al. 2014, 2015). The overlap be¬ 
tween the properties of Type II spicules and network jets has 
led to a growing conjecture that these represent identical mag¬ 
netic features that have been observed in different ways. It is 
clear that our modeled events with the highest speeds (corre¬ 
sponding also to the largest values of Az) have very similar 
properties as the observed features. 

Figure 5 highlights modeled events with Az > 2" (i.e., 
Az > 1.45 Mm) with larger and darker symbols. This di¬ 
viding line is close to the mean value of our distribution of 
simulated jet lengths, which exhibited a range between 0.1" 
and 7". The open BRAID flux tube has a diameter that ex¬ 
pands with increasing height, from about 1 " in the low chro¬ 
mosphere (zfe) to 3" at the largest heights shown in Figure 3. 
For comparison, the observed IRIS network jets have lengths 
between about 3" and 12", and roughly constant widths of or¬ 
der 0.5" (Tian et al. 2014, 2015). We believe that only events 
with Az greater than their cross-sectional diameter would ac- 



Figure 5. Two-dimensional diagram of representative time scales plotted 
versus vertical speeds. Observed ranges of parameters for Type I and II 
spicules (green solid-curve boxes) and IRIS network jets (red dotted-curve 
box) are compared with BRAID-model simulations of Af and Ejet. Large 
black symbols show modeled events with lengths greater than 2”, and 
small gray symbols show modeled events with lengths less than 2”. The 
Bertschinger & Chevalier (1985) relationship for periodic shocks is shown 
with a solid black curve. 

tually be observable as narrow jet-like enhancements in in¬ 
tensity. Thus, the darker symbols with Az > 2" are meant to 
show only the events that would be distinctly noticeable as 
jets in, e.g., IRIS image sequences. Shorter jets (like the mod¬ 
eled events shown in light gray) may exist on the Sun, but they 
would likely be buried in the rapidly fluctuating background 
of the underlying chromospheric network. 

The solid curve in Figure 5 shows the relationship between 
the initial upward velocities of shocks in a hydrostatic at¬ 
mosphere and the recurrence timescale between successive 
shocks in a periodic train. We solved the semi-analytic equa¬ 
tions given in Section II of Bertschinger & Chevalier (1985) 
for a range of Mach numbers between 1.1 and 30. In this 
model, gas is accelerated upward by each shock, then it tends 
to fall back down after the shock passes by. For periodic shock 
trains that obey the speed-timescale relation shown by the 
curve, each parcel of gas ends up at the same height one pe¬ 
riod later and executes cyclic motion. It is interesting that the 
properties of Type I spicules overlap with this critical curve, 
since they often appear to show parabolic trajectories in which 
parcels return to their original heights. 

Regions to the left of the critical curve in Figure 5 corre¬ 
spond to shorter timescales than are required for cyclic mo¬ 
tion. Thus, a given fluid parcel would encounter the “next” 
shock before it has fallen back to its original height. If the 
shock trains completely cover the stellar surface, such rapid 
recurrence would be associated with net mass loss (Willson 
& Hill 1979; Bertschinger & Chevalier 1985). Both Type II 
spicules and IRIS network jets occur in this region of param¬ 
eter space, and it has been suggested that these features feed 
plasma into the corona and solar wind (e.g., De Pontieu et al. 
2009; Tian et al. 2014). Chromospheric diagnostics of Type II 
spicules often show bright features moving up and not com- 
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ing back down (Pereira et al. 2012), but more recent coor¬ 
dinated observations with IRIS have shown some parabolic¬ 
like downflows (Skogsrud et al. 2015). Connections between 
observed features and the idealized periodic shock model of 
Bertschinger & Chevalier (1985) are instructive and sugges¬ 
tive, but they certainly do not tell the whole story. 

6. DISCUSSION AND CONCLUSIONS 

The goal of this paper was to explore one promising way 
that strong turbulence may produce dense, short-lived, and 
field-aligned extensions of chromospheric and TR plasma. 
We made use of a time-dependent RMHD model of Alfvenic 
turbulence in the coronal-hole network, and we took particu¬ 
lar notice of the intermittent amplitude variability in the mid¬ 
chromosphere. At a height of about 900 km above the pho¬ 
tosphere, the wave properties appear to be optimal to pro¬ 
duce a spike of nonlinear mode conversion into longitudinal, 
compressive fluctuations. These waves have been shown to 
be able to “puff up” the effective density scale height of the 
chromosphere and thus temporarily increase the height of the 
TR. Using an existing grid of models, we computed the time- 
dependent TR height as an instantaneous response to the vary¬ 
ing wave amplitudes and scale heights. Apparent upward ve¬ 
locities and recurrence timescales measured from the model 
time series agree quite well with the observed properties of 
IRIS networkjets and Type II spicules. 

There are other observable characteristics of jets and 
spicules that can, in principle, be compared with our mod¬ 
els. Do the IRIS jets appear only for specific ranges of plasma 
properties in the upper chromosphere (i.e., those that maxi¬ 
mize the mode conversion “spike” at Zb)l Is the observed Ail¬ 
ing factor of the jets (both in space and time) in agreement 
with the intermittency seen between the largest-amplitude 
pulses of (5v|| in Figure 4? These comparisons require robust 
statistics from the measurement of hundreds of individual jets 
and spicules. Collecting data with sufficient accuracy may re¬ 
quire the use of automated feature-detection algorithms (e.g., 
Aschwanden et al. 2013). 

We note that the models presented in this paper do not rep¬ 
resent a completely self-consistent simulation of the proposed 
jet/spicule formation mechanism. Instead, we attempted 
to show—via a sequence of separate, simple, and well- 
understood models—that the various ingredients are present 
at the right order of magnitude to produce the proposed ef¬ 
fects. We pay for this conceptual simplicity with the fact 
that the results (e.g., the up/down dynamics of the TR inter¬ 
face) are not likely to be quantitatively accurate. Full three- 
dimensional MHD simulations are required to test these ideas. 
Compressive simulations performed with fewer than three di¬ 
mensions (e.g., Matsumoto & Suzuki 2014; Kono et al. 2015) 
already show suggestive hints of the relevant mode conversion 
from Alfven waves to compressive/spicule-like pulses. 

A comprehensive explanation for the IRIS jets and Type II 
spicules will also require taking into account some additional 
processes and complexities that we did not include. For ex¬ 
ample: 

1. As described in Section 1, the larger polar plumes and 
jets are probably driven by magnetic reconnection at 
their footpoints. There may be an overlap between 
the smallest of these jets and the ones observed by 
IRIS in the (mostly) unipolar network (see also Moore 
et al. 2011). Coronal reconnection can also gener¬ 
ate Alfven waves (Hollweg 2006; Lynch et al. 2014), 


and turbulence can also lead to the formation of small 
scale reconnecting current sheets (e.g., Matthaeus et al. 
2015). The traditional dichotomy between wave-driven 
and reconnection-driven coronal heating theories is no 
longer so sharp, and it is useful to keep both kinds of 
processes in mind. 

2. Our models assumed that open magnetic flux tubes in 
the coronal-hole network are essentially isolated from 
one another. However, there has been a great deal of 
work to study how wave-like fluctuations can enable the 
sharing of energy between neighboring flux tubes and 
their weak-fleld surroundings (e.g., Uchida & Kaburaki 
1974; Roberts 2000; Bogdan et al. 2002; Hasan & van 
Ballegooijen 2008; Ofman 2009; Mumford et al. 2015). 
This is another type of “mode coupling” that needs to 
be taken into account. 

3. We ignored parametric instabilities and nonlinear mode 
coupling that is enhanced by inhomogeneous back¬ 
ground properties (see Section 2). The fact that we 
obtained the correct order-of-magnitude effect for the 
jet and Type II spicule properties may suggest that the 
adopted second-order ponderomotive coupling mech¬ 
anism is dominant. However, the other proposed ef¬ 
fects are likely to produce lower-frequency compressive 
waves that are the ones observed to survive to larger 
heights (see, e.g., Liu et al. 2015). 

Lastly, we note that the small jets and spicules discussed in 
this paper may be relevant to the larger problems of coronal 
heating and solar wind acceleration. There is evidence that 
the rapid upward mass transfer in Type II spicules continues 
as the plasma heats to temperatures in excess of 10^-10® K 
(McIntosh (fe De Pontieu 2009; Skogsrud et al. 2015). How¬ 
ever, there is also skepticism concerning the suggestion that 
jets and spicules act as a primary source of coronal plasma 
(Klimchuk 2012, 2015; Judge et al. 2012). Nevertheless, it 
appears more certain that the waves originating in lower atmo¬ 
spheric structures survive as they propagate up into the corona 
and heliosphere. Similar kinds of nonlinear mode coupling 
have been proposed to act along open held lines in the solar 
wind (e.g., Del Zannaet al. 2001; Chandran 2005; Cranmer & 
van Ballegooijen 2012; Miyamoto et al. 2014), and what we 
learn about this process in the chromosphere and TR can help 
improve our understanding of these other regions as well. 
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